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Abstract. We present a novel analytical approach for the calculation of the mean 
density of states in many-body systems made of confined indistinguishable and non- 
I interacting particles. Our method makes explicit the intrinsic geometry inherent in the 

symmetrization postulate and, in the spirit of the usual Weyl expansion for the smooth 
part of the density of states in single-particle confined systems, our results take the 
form of a sum over clusters of particles moving freely around manifolds in configuration 
space invariant under elements of the group of permutations. Being asymptotic, our 
approximation gives increasingly better results for large excitation energies and we 
formally confirm that it coincides with the celebrated Bethe estimate in the appropriate 
region. Moreover, our construction gives the correct high energy asymptotics expected 
from general considerations, and shows that the emergence of the fermionic ground 
state is actually a consequence of an extremely delicate large cancellation effect. 
Remarkably, our expansion in cluster zones is naturally incorporated for systems of 
interacting particles, opening the road to address the fundamental problem about 
the interplay between confinement and interactions in many-body systems of identical 
particles. 
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1. Introduction 

Understanding the general pliysical properties of interacting systems is one of the 
ultimate goals of classical and quantum mechanics, and probably the most difficult. 
An impressive toolbox including techniques such as the renormalization group [HE], 
perturbation expansions in Feynman diagrams [3|J1] and Density Functional Theory [SUB] 
has been developed during the past decades with the sole objective of constructing 
the energy spectrum of quantum systems where external potentials, inter-particle 
interactions and quantum statistics must be simultaneously considered. 

In a non-relativistic first-quantization scenario (which is an effective low-energy 
approximation to the fundamental field-theoretical relativistic description) both the 
Hilbert space and the Hamiltonian are fairly well understood and the problem of 
constructing the energy spectrum is reduced to the, still formidable, problem of 
calculating the many-body density of states 

g±{E, N) = [tri^^ G{E + ie)] . (1) 

TT L J 

Here is the total number of particles, (±) refers to the symmetry with respect 
to particle exchange of the subspace where the trace operation is performed (fully 
symmetric wavefunction for bosons and fully antisymmetric one for fermions), and G 
is the Green function of the many-body system, defined in terms of the full interacting 
Hamiltonian 

H = Ho + gU (2) 

as ^ 

G{z) = {z-h)~ , (3) 

where Hq describes a set of non-interacting particles, U is the (two-body) interaction 
potential and g the coupling strength. 

In this paper we are interested in the situation where the quantum system is made 
up of well defined microscopic degrees of freedom either bosonic or fermionic, subject 
to a hard wall confinement in D dimensions such that the system is bounded for all 
energies. In this case, the density of states has the generic form 

Q^{E,N) = Y^6{E-E\^^), (4) 

n 

where n = 1,2... labels the ordered eigenvalues e[^\E2^\... of H for fixed total 
number of particles A^. Although the precise knowledge of the density of states provides 
all the information about the spectrum of the system, general considerations show that 
g±{E, N) can be unambiguously decomposed into a smooth and an oscillatory part, 

g±{E, N) = g±{E, N) + g"^'{E, N) , (5) 

where the scaling of the smooth part is given asymptotically by the Weyl formula 

^±(E,Ar)~l/n^^. (6) 
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If the interaction and the confinement is such that no constants of motion remain besides 
the total energy and the classical many-body dynamics is chaotic, for the oscillatory part 
we can formally apply the Gutzwiller trace formula (IDlIIl] to get 



in the semiclassical limit ^ — )■ Ij. From the scaling with h in the semiclassical limit 
it is expected that oscillatory contributions to the many-body density of states in the 
strongly interacting regime where the quantum mechanical description of the system lies 
beyond the single-particle picture, are extremely small compared to the mean density 
Q±{E, N). Understandably, continuous effort has been dedicated to develop methods to 
construct this function either numerically or analytically [H[TT |[T3]4T5] . 

To date, the only general way to obtain the precise g±{E,N) is by explicit 
(typically numerical) diagonalization techniques of the many-body problem followed 
by a convolution with a smoothing function We{E) of width e. 



However, due to the fast (exponential) growth of the basis required to achieve good 
convergence of the numerical results, such numerical approaches can deal with only 
moderate numbers of particles (see, e.g., |16] for state of the art calculations). 

Given the complexity of the problem, alternative methods are mandatory. Self- 
consistent mean- field methods, firmly grounded in the Kohn-Sham theorem [HE], 
provide the most efficient way to construct a set of single-particle wavefunctions 
such that the ground state energy of the interacting system can be systematically 
approximated by artificial single-particle energies supplemented with the appropriate 
symmetry of the many-body wavefunction. 

Despite the extremely successful application of self-consistent methods, ranging 
from nuclear physics to molecular systems, reaching chemical problems [HlEllIT], the 
calculation of the many-body density of states within mean-field approaches faces a 
conceptually deep and basically unsolved problem. It is related with its very definition, 
and stems from the fact that the calculation of different types of observables actually 
requires a different definition for the mean field. Examples of this ambiguity are the 
calculation of ground-state vs excited-state energies and the construction of static vs 
dynamical properties of the system (for the calculation of transition amplitudes, for 
example, the mean field is necessarily time- dependent and depends on the initial and 
final states as well) [3]. 

To be more precise, extending self-consistent approaches like the Density Functional 
Theory or Hartree-Fock method to excited states above the ground state requires 
a specific knowledge of which single-particle orbitals are optimized through the self- 
consistent equations, and therefore each excited many-body state requires a separated 

I We do not discuss questions [12] concerning the derivation of ([7]) for DN > 3 here. 



(7) 





A geometrical approach to the mean density of states in many-body quantum systems 4 

calculation leading to a different mean-field. Similar problems appear in formulations of 
the many-body problem based on a functional representation of the propagator where 
the mean-field is defined by means of a saddle-point equation: roughly speaking, the 
single-particle artificial potential used to mimic the effects of interparticle interactions 
becomes dependent on the excitation energy itself [ij. 

Working with a mean-field which is simply a function of the position, independent 
on both excitation energy and/or time is then a strong assumption that lies behind much 
of the efforts to understand many-body systems in terms of a single-particle picture, 
an assumption which is in fact turned into essential in most cases where even a mild 
dependence of the mean field on the excitation energy would make the calculations 
impossible. 

Luckily, many important physical effects accessible to experimental observation take 
place near or at the ground state energy, and it is expected that in this situation the 
mean field, providing the independent particle picture for the ground state within any 
of the self-consistent methods, can be used to calculate physical properties of excited 
states with decreasing accuracy as we move to higher and higher excitation energies. 
This physical consideration, together with the more pragmatic reasons explained above, 
has shown to be remarkably useful. In fact, the idea of an unique mean field allows us 
to use all the standard machinery of single-particle physics as input for the statistical 
mechanics results for independent particles and to finally produce experimentally 
accessible predictions beyond the strictly non-interacting case. A good mean field is 
then an excellent starting point. 

A paradigm of the success of this approach is the study of low energy excitations 
of bounded fermionic systems with many particles like nuclei, metallic clusters and 
quantum dots [5],[6lll7j. Here, a self-consistent calculation of variational type is set up 
to fix once and for all a single-particle potential responsible to replace the interparticle 
interaction. The result of the numerical calculation is used then to fit the parameters of a 
large family of functional forms (e.g. Woods-Saxon for nuclei [T8H20]) to finally produce 
an analytical form. Once this is achieved, standard methods of statistical physics are 
applied. 

A result of these studies which is of importance for the present work is that, for 
a large enough number of electrons interacting through Coulomb forces in a billiard- 
shaped quantum dot, the mean field (in the sense of Hartree-Fock) closely follows the 
billiard potential. This allows us to assume a billiard model for the interacting system 
in mean field approximation knowing that it already captures part of the physics beyond 
the non-interacting case through a softening of the boundaries. 

The program outlined above has reached a high sophistication, in particular when 
the single-particle physics is treated analytically by means of semiclassical methods, 
well suited to study the effective single-particle problem around the Fermi energy for 
the regime of large numbers of particles. The idea here is to construct the single-particle 
spectrum in terms of periodic orbits of the classical system considered as a single-particle 
moving within the mean field potential. The amount of work on the subject is huge and 
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we refer to [H] for a more complete exposition including experimentally relevant effects 
like the existence of shell effects visible in the many-body density of states and due to 
continuous families of classical periodic orbits. 

Another possibility offered by the mean field approach comes from one of the main 
lessons we have learnt from single-particle systems, namely that contrary to g±{E,N) 
which would be sensitive to every single detail of the microscopic mean field Hamiltonian, 
the smooth part of the density of states g±{E,N) predominantly depends on few 
classical quantities related to the measure of the classical phase space manifolds at 
given energy. Therefore, instead of going through numerically expensive calculations 
in order to construct the single-particle energies in mean field approximation just to 
average out again the density of states and get its smooth part as indicated by ([8]), one 
can try to construct and understand the classical phase-space structures behind . This 
approach can be carried on in two different and complementary ways. 

The first option is to push the mean field picture and calculate the smooth part 
of the single-particle density of states using the Thomas-Fermi approximation and its 
variants [11] , and use it together with well estabhshed statistical techniques to construct 
g±{E, N). This is the line of thought that lies behind the original attempt of Bethe [13] 
to calculate the smooth part of the level density for nuclei and, more recently, a similar 
concept is followed in the seminal works of Weidenmiiller [21] in order to formally 
construct the exact level density g±{E,N) for non-interacting particles. It is important 
to mention that this approach is based on the classical single-particle phase space, and 
the construction of the many-body density of states is purely formal and has no direct 
interpretation in terms of the classical many-body phase space. 

The second possibility, namely, to calculate the smooth part of the many-body 
density of states in mean field approximation by relating it directly to the structure of 
the many-body classical phase space has not been systematically addressed before. In 
our opinion, such approach has the potential advantage of avoiding the a priori conflict 
with the inclusion of residual interactions inherent to any approach that is based on the 
single-particle phase space. 

As we will show, these two approaches give quite different points of view for the 
calculation of g±{E,N), and their equivalence for the strict mean field limit is by no 
means trivial. However, once this equivalence is established, our methods based on the 
many-body phase space will allow to address the fundamental question concerning how 
the relevant geometrical structures are translated into the problem in the presence of 
residual interactions. This is the program we propose here. 

Once the single-particle picture is adopted the study of reference is the so-called 
Bethe estimate, providing an asymptotic result for the density of states in many-fermion 
systems in mean field approximation, valid for energies far enough (in units of the 
single-particle mean level spacing) from the ground state Eqs and large numbers of 
particles [13], 
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The essential aspect of Bethe's result in the context of interest here is that it is 
an asymptotic approximation for the smooth part of the density of states, and can 
be interpreted as the fermionic analogue of the Weyl expansion for billiard systems. 
However, and despite its enormous importance, (Q (and the thermodynamic formalism 
used in its derivation) is of limited use for us, as it does not provide any clue on 
how classical phase space manifolds work together to produce both its characteristic 
functional form and the scale -Eqs- 

In the physical scenario of interest here (the application of Bethe's method in 
confined electronic systems) the most relevant problems besides the all long issue of 
residual interactions are, i) the inclusion of finite N effects, ii) the consistent treatment 
of the ground state energy and the density of states around it, and iii) the emergence 
of the standard Weyl expansion q±{E,N) ~ E^^^"^"^ for high enough energies. These 
questions have been addressed already, and the amount of literature in the subject is 
extensive, so we will only briefly review the state of the art. 

Finite size effects on the asymptotic density of states in systems of fermions can be 
systematically calculated in mean field approximation by extending Bethe's result, which 
is the leading order in an expansion valid for large excitation energies and large numbers 
of particles obtained by a saddle point approximation of the exact grand-canonical 
partition function [^ [T3HT5] . Because the problem of counting many-body eigenstates 
for non-interacting identical particles has a natural combinatorial formulation, this 
approach has produced an unexpected and fruitful interaction with number theory. 
In this spirit, corrections to the Bethe estimate arising from finite number of particles, 
oscillatory corrections to the single-particle density of states, and shell effects affecting 
the ground state energy have been considered [22] . 

The status of the ground state energy within the asymptotic approach is somehow 
delicate, as strictly speaking, the many-body density of states must vanish identically 
for energies bellow Eqs- As obvious as this observation may be, it turns out that it 
is extremely difficult to construct a theory providing the large energy asymptotics for 
g±{E,N) while keeping this condition exact. This problem may be considered at first 
glance a merely academic, since after all no physical process takes place for energies 
bellow the ground state, and the later has a very precise definition in terms of the 
single-particle density of states. We must however keep in mind that the final goal of 
any approach is to deal with the fully interacting system and/or to provide insight and 
better methods to define and calculate the mean field and the effect of the residual 
interactions. Beyond the mean field picture, we eventually need a systematic method to 
identify and construct the ground state energy independent of the counting prescription 
valid in the non-interacting case. To the best of our knowledge, a systematic study of 
how do approximations for g±{E, N) behave for E < Eqs is missing. 

Finally, very general and robust considerations demand that asymptotically (when 
the energy goes to infinity), essential quantum mechanical effects such as the non-zero 
ground state energy gradually disappear and a purely classical description emerges. In 
this limit, one should recover the standard Weyl expansion for the density of states 
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where quantum symmetry effects only appear as a global reduction of the available 
phase space to a fraction of due to the identity of the particles but not to their 
particular statistics [23j . Using asymptotic methods to understand this transition leads 
to some interesting results connected with the number-theoretical formulation of the 
problem [22]. 

Together with the goal of providing a geometrical approach to g±{E,N), the last 
paragraphs indicate the main motivations of the present work. We attempt to provide 
a method to construct the Weyl approximation to the smooth part of the density of 
states which relies only on kinematic and geometrical aspects within the mean field 
approximation. As expected, our results are connected with several others (in particular 
with Bethe's) and an important aspect of our work is to make these connections explicit. 
However, the method itself and the physical idea behind it, namely that the Weyl 
expansion can be systematically constructed out of free propagation near symmetry 
manifolds are our novel contributions to the subject. 

The paper is organised as follows. In section [2] we introduce the basic notation and 
briefly review the construction of the Weyl expansion for systems of non-interacting 
identical but distinguishable particles. In sections [3] and H] the symmetrization postulate 
is used to give the formal expression for the full density of states for indistinguishable 
particles, and we show how this formal object can be understood in terms of the 
geometry of a higher dimensional phase space when the fundamental domain associated 
with the group of permutations Sn is considered. The role of classical manifolds 
invariant under different elements of the symmetric group is clearly seen using the 
example of two fermions on a line in section |5l The very non-trivial generalization of this 
construction for arbitrary type of particles (bosons or fermions), arbitrary dimension D 
of the single-particle configuration space Q and arbitrary number of particles is fully 
carried out in sections [6] and [7] and culminate with a full identification of the classical 
manifolds and their measures responsible for the functional form of g±{E,N) and the 
emergence of Eqs- We analyse our results in sections [8] and |9l where the equivalence 
of our results with the Weidenmiiller convolution formula and with Bethe's estimate 
are rigorously proved. We conclude and discuss the extension of our results for the 
interacting case in section [101 

2. Non-interacting Distinguishable Particles 

First consider a billiard system of non-interacting, identical, but distinguishable 
particles in D dimensions, specified by their coordinates 



D 



I = 



l,...,Ar 



(10) 



with spatial components 




d= 1, ...,£). 



(11) 
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This corresponds to an effective ■ Z)-dimensional billiard system of a single particle 
described by 

q=(qi,...,q;v)GM^^. (12) 

In this case the only effect of particle exchange symmetry is thereby the existence of 
discrete spatial symmetries in the effective higher dimensional single-particle system. 
But due to distinguishability, there is no restriction to any subset of wavefunctions with 
specific symmetry under exchange transformation. 

In general, the density of states (DOS) of a bound, time independent system 
can be expressed as the inverse Laplace transform of the trace of the propagator 
U{t) = exp{-i/hHty. 



e{E) = 



J d^^q K{q, q, t = -ih/3) 



(E) , (13) 



where the trace is performed in coordinate space with K{q', q; t) = (q'| f/ (t) |q) and the 
inverse Laplace transform has to be applied with respect to the variable /3. 

In semiclassical approximation (corresponding to the formal limit h ^ 0) there are 
two contributions to the DOS 

g^^\E) = g^^^{E) + g{E), (14) 

one oscillating {g°^'^) and one smooth (g) in the energy E. The oscillatory part arises 
from various stationary phase approximations in (fT3!) starting from a path integral 
representation of the propagator. The process leads to a description by periodic orbits of 
the underlying classical system. The oscillatory part of the DOS is then in semiclassical 
approximation expressed by the Gutzwiller trace formula in the chaotic case p^[2^ and 
the Berry- Tabor trace formula in the integrable case [25] respectively. 

The smooth part of the DOS is related to short path contributions that are not 
caught by periodic orbits in the analysis of the trace of the propagator. Reflecting the 
short time behaviour of the propagator, these are related to the assumption of local free 
quantum propagation and additional boundary corrections [26]. The Weyl expansion 
for the N ■ D-dimensional billiard reads 

m \— Vnd 



^ 2 ^ (15) 



-. ND-l n 



The first term in f|T5l) will be referred to as the volume Weyl term and equals the 
Thomas- Fermi approximation ^tf(-E') [H] proportional to the A^-D- dimensional volume 
Vjvu- The second term originates from wave reflection near billiard boundaries under 
the assumption of local flatness of the surface. This involves free quantum propagation 
to mirror points yielding a fast converging integral over the coordinate perpendicular to 
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the surface in ( I13p . Thus the second term is proportional to the integral over parallel 
coordinates yielding the surface Snd-i (not to be confused with the symmetric group 
Sn) of the billiard instead of its volume. Higher corrections in the expansion correspond 
to propagation between multiply reflected image points accounting for curvature, edges 
or corners of the boundary. 

3. Non-interacting Undistinguishable Particles 

In the case of identical particles that are undistinguishable, the state of the system 
obeys a specific symmetry with respect to particle exchange. The state is either 
symmetric or antisymmetric under the exchange of any two particles, depending on 
whether they are bosons or fermions. 

P|^±) = (±1)^ , (-1)^ := sgn(P) (16) 

for any permutation P G Sn, where P is the corresponding permutation operator acting 
on many-body states, and plus and minus refer to bosons respectively fermions. In order 
to obtain the physical spectrum of such a system one has to restrict to those eigenenergies 
that correspond to the subspace of Hilbert space with appropriate symmetry. Let 
1± = ij_ be the projector onto the subspace of correct symmetry. Restricting the 
trace in f[T^ to this subspace is equivalent to replacing the propagator by its symmetry 
projected analogue 

U±{t):=i±U{t)U = l±U{t), (17) 
^±(q', q; ^) ^= M E (±1)'^ ^(^q'' t) ■ (18) 

■ P^Sn 

In f|T7|) . the commutation of the time evolution operator and the symmetry projector 
due to [P, H] = and idempotence of 1± have been used. This leads to the symmetry 
projected DOS 

Thus symmetry causes the need of taking wave propagation over finite distance into 
account, as Pq 7^ q in general. For the oscillating part in g±{E) it is possible to construct 
a fundamental domain in phase space where it is again sufficient to find periodic orbits 
and additionally the group characters (±1)^ of the group elements P connected to each 
trajectory (q, p) 1— )■ (Pq, Pp) in the unfolded full domain [27|. However this procedure 
has no direct application to the smooth part of the DOS in the general case of arbitrary 
particle number and spatial dimension. 

This is true even in a bosonic system in D > 1 despite the possibility of 
mapping it to the higher dimensional system of a single particle without symmetry 
moving in the fundamental domain in coordinate space with topological identification 



1 

m 



EM 

Pe5w 



d^^qK{Pq,q;t = -ih/3) 



(19) 
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of symmetry related points. That is because of the non-trivial structure of such a 
wrapped fundamental domain especially in the vicinity of symmetry planes defined by 
Pq = q for some P G Sn- In order to illustrate this, compare to a two dimensional 
single-particle system with discrete rotational symmetry 

R^ = exp(^-^L,^^ , (20) 

= — , n G N. 

n 

The restriction to the subspace of states symmetric under the elementary rotation 
R(l> iV'sym) = l^sym) IS equivalent to the restriction to the wrapped fundamental domain 
with usual wave dynamics. The wrapped fundamental domain in this example is the 
restriction of the billiard to a sector of central angle with identification of points along 
its two bordering half-lines. As it is equal to a cone, this produces non-trivial wave 
propagation at the origin which gives rise to additional corrections in the level density 
QsymiE) of symmetric states [28]. Analogue to that, mapping a bosonic system to its 
wrapped fundamental domain implies non-trivial wave propagation in the vicinity of 
the symmetry planes. This shows that it is reasonable to stay in the full domain for the 
calculation of the smooth part g±{E), and this is the approach we will follow here. 

Previous to the treatment of the general case it is instructive to analyse the simple 
example of many identical fermions on a line which can be mapped to a fundamental 
domain where the additional correction due to symmetry can easily be obtained by usual 
methods. 



4. Equivalence of Many-Body and High Dimensional Single-Particle 
Pictures 

This section will mainly focus on systems of many fermions moving on a line of length 
L. These systems have some special properties that are setting them apart from higher 
dimensional ones. It is worth restricting to such systems for a moment since they can 
easily be mapped onto single-particle systems. 

In one dimensional systems a fundamental domain in coordinate space can be given 

by 

^ := |q G < gi < g2 < ■ ■■ < gTV < i^} • (21) 
Its boundaries are given by the equations 

= qi+i , z = 1, . . . , iV - 1 (22) 
due to symmetry related reduction and 

gi = 0, 

^ ' (23) 
qN = L 
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Figure 1: Fundamental domain ^ for N = 3, D = 1, bounded by symmetry planes 
qi = q2 and q2 = qs (light grey) and physical boundaries qi = and q^ = L (dark grey 
and shaded). 



for the physical confinement to the line [0,L]. An example of this construction in the 
case of three fermions is shown in figure [1] The full domain can be reobtained by 
applying all possible permutations P G Sj\f to the fundamental domain 

[0,L]^= [j P{^). (24) 

Due to D = 1, topological identification of boundary points in this context is not needed. 

For fermions the restriction to antisymmetric states yields the condition of vanishing 
wave function all along the boundary 

^(q) = (Vq)(3z^j)(?. = gj)- (25) 

As a Dirichlet boundary condition, this condition is sufficient to determine the 
eigenf unctions in ^ together with the single-particle conditions f l23|) . The symmetry 
planes ( 122|) can be thought of as hard walls or in other words infinite potential barriers. 
The values of the wave functions in the other parts of the full domain are then obtained 
by 

V'(Pq) = (-l)^^(q) . (26) 

Thus a ID billiard with fermions is equivalent to a single-particle billiard of dimension 
N ■ D, in which the usual Weyl expansion can be used to obtain g^{E). 

One has to stress that this is a feature of one dimensional systems only because 
there no additional condition besides ( l25l) is imposed on the wave function within the 
fundamental domain. In contrast to that, the corresponding condition ip{q) = for 
dimensions larger than one is not given along the whole boundaries of the fundamental 
domain, but instead only on lower dimensional manifolds embedded in those. The reason 
is that for D > 1, ^ is not bounded by the symmetry planes q^ = q^, which have a 
dimension oiN-D — D<N-D — 1 and therefore are not able to separate volumes 
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in N ■ D-dimensional configuration space. Rather the restriction to one of the spatial 
components in the conditions that define the symmetry planes are able to do so. One 
choice of fundamental domain is 

^:={q6^]^|g«<g?)<•■■<g«}, (27) 

where its symmetry related boundaries are defined through 

g«=gJi\, z = l,...,iV-l, (28) 

and q G fi^ denotes the physical confinement to the interior Q C of the billiard. 
The symmetry planes = q^+i, on which ( l26ll imposes vanishing of the wave function, 
are only lower dimensional submanifolds embedded in the boundaries. Thus, the 
condition (1261) can not be expressed as a Dirichlet boundary condition in a fundamental 
domain. 

Moreover, even in the case of D = 1 sharp edges in the boundary of ^ can cause 
problems for N > 2, making the expansion of Balian and Bloch [26] inapplicable, which is 
a generalization of the Weyl expansion to arbitrary dimensions but smooth boundaries|§| 



5. Two Non-Interacting Fermions on a Line 

Consider a one-dimensional system of two fermions confined to a line of length L. The 
only two permutations are the identity and their exchange. Figure [2a] illustrates the 
two possible contributions to the propagator. Since the system has an effective two- 
dimensional description it is straightforward to compare it to a single-particle two- 
dimensional billiard. There, g{E) is made up of contributions from free propagation 
and from reflections on the boundary, which is illustrated in figure [20 

As discussed in section [H there is a simple two-dimensional single particle billiard 
that is exactly equivalent to the two particle system. That is, the billiard defined by 
the fundamental domain, here chosen as ^ : L > gi > g2 > with an additional 
hard wall boundary along the symmetry line qi = q2- In this two-dimensional picture 
reflections on the additional boundary, which are addressed by propagation to mirror 
points, are mapped to the propagation with respect to the exchange permutation in the 
one dimensional many-body picture. Figure [3] illustrates the corresponding propagation 
in both pictures. The smooth part of the DOS of the two-fermion system including 
corner corrections reads 

*-(^) 4 if) "f^) - p + K if) ' ^-*''(^' + >^ • 

§ It is worth to note that one could topologically identify points along the boundary of ^ that are 
related by symmetry and thereby create a fundamental domain with complex topology. We will call this 
object the wrapped fundamental domain. The problem in the fermionic case then would be the loss of 
continuity of the wavefunction because of sign inversion ijj (q) —rp (q) in the direction perpendicular 
to boundaries related to odd permutations. This condition seems quite peculiar and so far the author 
has not found a treatment of it in [26] . 
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(a) 2 particles in 1 D (b) 1 particle in 2 D 



Figure 2: Comparison of the system of two identical particles in one dimension and one 
particle in two dimensions. 




(a) 2 particles in 1 D (b) 1 particle in 2 D 



Figure 3: Equivalence of the system of two identical particles in one dimension and one 
particle in the two-dimensional fundamental domain. 

The Weyl expansion fl29p includes a volume term Qv{E) with area A = L'^/2, where the 
factor of I originating in the restriction to the fundamental domain corresponds to the 
factor in the symmetry projected propagator ( IT8l) . Thus this factor in the volume 
term, which corresponds to taking into account only the identity permutation, is the 
leading order effect of exchange symmetry. In general exchange corrections related to all 
other permutations are of sub-leading order as they correspond to Weyl-like boundary 
corrections. For example, in the case at hand, the only exchange permutation yields 
in leading order the perimeter correction proportional to a/2. Nevertheless, in general 
all exchange corrections are important to give physically reasonable results, as will be 
shown in the next sections. 
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6. General Case - Propagation in Cluster Zones 
6.1. Invariant Manifolds and Cluster Zones 

The last section showed that in the example of two particles on a line the correction to 
Q-{E) due to the exchange permutation is related to the propagation in the vicinity 
of the symmetry line qi = q2 just as the inclusion of wave reflections in a single- 
particle billiard only affects the short time propagation near the physical boundary. The 
symmetry line is characterised by the invariance under P = (12), so that the distance 
|Pq — q| becomes zero. This is the very reason to assume short path contributions to 
come from its vicinity. The concept of invariant manifolds is extended to the general 
case by finding the manifolds associated with each permutation P, defined by 

Mp = {qeR^^ \\Pq-q\=0}nQ^ . (30) 

P can be written as a composition of commuting cycles (see for example [29] ) 

P = ai---ai, (31) 

acting on distinct sets of particle indexes of size 

N,,N2,...,Ni. (32) 

So we see that A4p is the manifold defined by the coincidence of the coordinates of all 
particles associated with each cycle 

I 

7Wp= fllqefi^ I q. = q,- Vz, j G . (33) 

UJ = 1 

As a simple example take the permutation 

P=( 134)(25), (34) 

whose associated manifold A4 p corresponds to the condition (see figure Hj) 

(qi = qs = q4) a (qa = qs) . (35) 

In contrast to the symmetry line in the one-dimensional two particle case, these 
manifolds in general can not be seen as a boundary or surface in coordinate space 
in the sense of dividing the space into distinct pieces. The vicinities of these invariant 
manifolds shall from now on be referred to as cluster zones. All particles associated with 
a particular cycle index subset will be subsumed to the notion of a cluster. A system 
that is momentarily arranged in a particular cluster zone is composed of / clusters, each 
associated with a cycle in P. Each cluster u is composed of A''^^ particles according to 
the length of the cycle. 
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1^4 ^2* ^ ^q2,5 

F=(134)(25) Mp 

Figure 4: Example of the correspondence between the cycle decomposition of a 
permutation P and the invariant manifold Aip pictured as the associated clustering 
of particles. 



A separation into coordinates parallel to A^p and perpendicular to it suggests itself, 
since the propagation near M.p will not depend on shifting the position along them. 
This holds at least as long as one does not get too close to the single-particle bilhard 
boundaries to be included later. First, these single-particle boundary corrections shall 
be neglected, assuming the propagation to be invariant along M.p. This will be referred 
to as the unconfined case. 

Furthermore, the invariance of propagation along the invariant manifolds is in a 
strict treatment also broken in the case of interaction. But the intention lying behind 
this construction is that when restricting to interactions of rather short range, the 
propagation can be assumed to be invariant along A^p as long as one does not get 
too close to other invariant submanifolds. In other words, as long as the coordinates 
corresponding to different cycles do not become too close. Or to put it into a more 
intuitive picture, the different clusters should not collide. A discussion on that can be 
found in the last section. 



6.2. The Measure of Invariant Manifolds 

The infinitesimal volume element d/i on A^p is determined by the infinitesimal vectors in 
full (A^D)- dimensional coordinate space lying in A^p that correspond to the variations 
of independent coordinates. 

First, for each permutation P in cycle decomposition ( l3TI - [32l) the particles are 
relabelled without loss of generality in a way that the first cycle ai involves the first A^i 
particles, the second cycle (J2 involves the A^2 subsequent particles, and so on. 

Since all particles associated with one cycle have to fulfil the condition of equal 
coordinates f p3|) . there is exactly one independent Z)-dimensional vector for each cycle 
a^j , OJ = I, . . . ,1 , while all other particles in the same cycle have to follow in order to 
remain on Alp. Let the / independent vectors be denoted by 

x^ = (x«,...,x£^))Gr], u = l,...,l. (36) 

With these preliminaries, any point on Alp is described by 

q = (xi, . . . , xi, X2, . . . , X2, . . . , X,, . . . , X/) e . (37) 

JVi N2 Ni 
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The infinitesimal volume element d/i of A^p is the volume of the parallelotope spanned 
by all I ■ D infinitesimal tangent vectors 

dx^f , w = 1, . . . , /, d = l,...,D. (38) 



dx, 

(dyLi)^ equals the Gramian determinant G of all these vectors, which is the determinant 
of the matrix made up of all pairwise scalar products. From (J371) one obtains 



(9q 

dx^J^ 

with 



: 0,...^..,0 ,e,,...,erf, 0,......,0 ) (39) 



erf = (0,... ,0,1,0,. ..,0) gM^ (40) 

t 

d-th 



and therefore 

dq dq 



^ dx[f' dx^^, ■* 
The Gramian determinant reads 



Sdd'Sujuj'N^ ■ (41) 



' 2 



G' = det(diag(iVi,...,iVi,iV2,...,iV2,...,iV;,...,iV0) ■ (jj d^'x^) > (42) 

D D D "^=1 



dfi = VG= v^A^ ■ ■ ■ ■ d^xi ■ ■ ■ d^xi . (43) 

The total measure of the manifold fi{Aip) is obtained by integration of all independent 
coordinates x^^ over fl: 

^^{Mp) = fdfi = Vi,- (j[N^y , (44) 

with the dimensional volume of the billiard 

Vd = J d^x^ . (45) 

This measure is only depending on the partition of into integers A^i + ■ ■ ■ + A*"; 
corresponding to the decomposition of P into cycles with lengths Ni, . . . , Ni. Each 
permutation associated with the particular partition {A^^i, . . . , A*";} yields an invariant 
manifold of the same measure (HI]). Furthermore the contributions from short path 
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propagation in their vicinities will be the same due to the symmetry with respect to 
relabelling particle indexes. This allows the replacement of the sum over all permutations 
by a sum over all distinct partitions of with an additional factor of 

I 

c(N.,...,N,) = m[Xlj^){^ ) (46) 

w=l n=l 

in each summand, which is the number of permutations P G Sn with cycle lengths 
{iVi, . . . , iV/} in their cycle decomposition. Thereby, m„ denotes the multiplicity the 
cycle length n appears with. 

7. The Weyl Expansion of Non-Interacting Particles 

Up to this point, the analysis is quite general and is valid also for interacting systems. 
Although we believe it to be feasible in the context of interaction, for now the non- 
interacting case will be carried out explicitly. Furthermore, in this calculation, effects 
of the physical boundary are omitted. 

Under these assumptions the propagator in (fTSjl is taken as the product of free 
propagators of all particles 

N 

ir(q',q;t)= J]iro(qU^;^)• (47) 

i=l 

At the end of this section, the full expression including effects of (locally flat) physical 
boundary can be found. The corresponding calculation is shown in the appendix. 

For the calculation of the summand corresponding to a particular permutation P 
in (fT9l) . the particle indexes of different cycles in P don not mix up, yielding a product 
of independent propagators, which can be traced separately, each factor corresponding 
to a specific cycle in the decomposition ( !3T1) . Consider now the trace of all coordinates 
corresponding to the cycle a^. For this purpose, the particles are relabelled, so that the 
indexes associated with a^^ are simply /^^ = {!,... N^}. For the sake of simplicity we 
write q and Pq by meaning the restrictions to the first particles 

q= (qi,q2,...,qn), ^^^^ 

Pq = (qa, . . . , qn, qi) • 

Furthermore note the abbreviation n = N^. The integral over the associated coordinates 
reads 

//I TlD 
d-,/.„(Pq.q;<) ^ / d"", (^) " exp(l - q^) . (49) 

where the equality of a product of n free propagators in D dimensions and one n ■ D- 
dimensional free propagator has been used. With the distance vector 



^q - q = (q2 - qi, qs - q2, • • • , qi - qn) , 



(50) 
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the squared distance is 

|Pq - = |q2 - qiT H h |qi - q^l^ 

=f:[i<i?-iy+---+i9r'-m- (51) 
d=i 

The overall squared distance is the sum of squared distances according to one spatial 
component, which are just the summands in (I5T1) . The following calculation proceeds 
in equal manner for all spatial components d = 1, . . . , D. So the notation is further 
simplified by calculating only the factor corresponding to one spatial component in ( H9l) 
and by omitting the superscript (d). For this calculation in n-dimensional space we 
simply write q and Pq by meaning the corresponding tuples of one particular spatial 
component. 



-Pq= (92,g3,---,gn,gi). 

Which of the definitions f lT2|) . f l48|) and (|52|) is used in a particular step will be clear 
from the number of regarded dimensions in the context. 

In this simplified notation (152|) each summand in (I5T]) is 

|Pq - q|2 = - + . . . + (g„ - g„_i)2 + (g^ _ g^f . (53) 

The trace to calculate is 

/ d".A.(^0.<.;*)^/d". {^f exp(i^lPq-qP). (54) 

The squared distance is of second order in all coordinates, enabling to perform the 
integral fl54p as a generalized multidimensional Gaussian integral 



I d™x exp (-^x^Mx) = ^ , M = M^eGL^, (55) 

which will not be used directly, since the determinant of M equals zero. It has one 
eigenvalue A = corresponding to the direction parallel to the invariant manifold, 
expressing the local translational invariance along 

q|| = ^(l, •••,!), (56) 



since the distance vector is invariant in this direction, 

P(q + aqy) - (q + aq|[) = Pq - q + a(Pq|| - q||) = -Pq - q • (57) 
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This is the very reason why the separation into coordinates parallel and perpendicular 
to the invariant manifolds suggests itself. One way to proceed would be to introduce 
suitable perpendicular coordinates, perform the corresponding lower dimensional 
integral of the propagator and in the end multiply it by the measure of the invariant 
manifold. This would be the direct analogue to the usual computation of the surface 
correction in the single-particle Weyl expansion f|T5l) . One has to stress that in the 
interacting case this is most likely the most convenient way to calculate the trace of the 
propagator. And indeed, one can follow this procedure in the non-interacting case. But 
the introduction of perpendicular coordinates is rather uncomfortable since naturally one 
is lead to non-orthogonal coordinate systems and therefore has to introduce a metric 
tensor and pay extra attention to the arising volume elements. Although the calculation 
for the free case can be carried out in this manner, we will present an alternative 
approach in the non-interacting case that is more convenient and straightforward. 

For the following analysis, a minimum cycle length of n > 2 is assumed. The 
trivial case n = 1 will be included automatically in the resulting expressions. In the 
n-dimensional space of particle coordinates corresponding to only one cycle and only one 
spatial dimension, the subspace of vectors under which the squared distance is invariant 
is only one-dimensional (there is only one qy). Accordingly, the matrix M has exactly 
one eigenvalue that is vanishing when bringing the trace (l5ll) into the form of (155|) . 
Therefore it is sufficient to separate one of the n coordinates, e.g. qi and calculate the 
integral over all others as a generalized multidimensional Gaussian integral with linear 
term 



d". exp(-ix-Mx + v-x) = exp(iv-M-vj (58) 

with V G C™" and M = M'^ G GL„i. The remaining integral Jdgi can then 
be kept and eventually, when considering all spatial components and cycles, it will 
automatically produce the measure of A^p together with the determinant prefactors. 
Parts of the prefactors will thereby act as the Jacobian determinant associated with 
the relation of the volume element of the manifold to the independent coordinates q[^l 
, d = 1, . . . , D ,u = 1, . . . ,1. We abbreviate 

and write (|5^ as 

(~^)' / dgiexp(2agi) / dga ■ ■ ■ dg„ exp(^-^ ^ Aij-gi+ig^+i a ^ (60) 

ij=l i=l 



with some symmetric matrix A and a vector b, which are identified by separating all 
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gi-dependent terms in (1351) : 



n-l 



|Pq - q|2 = 2ql + ^ 2q^ - ^ 2g,g,+i - 2q,q2 - 2q 



i=2 



i=2 



A, = -4 



hi = -2qi{6ii + , 
A is a tridiagonal matrix of dimension n — 1 



z,j = l,...,n- 1. 



/ 1 -| 



A = (-4) ■ 







1 



-i 



whose inverse and determinant are 



i < J, 



det{A) =n(-2)"-^ 
Using ( ETh and (1^ gives 

= -Aql. 

With this and ( 158|) the whole integral (!60|) becomes 



(27r)"-i 
Try Y a" -Met (A) 



dgiexp(2ag^H — aYy'' A ^b 



(61) 



(62) 



(63) 
(64) 



(65) 



— — j ^ n 2 I dgi . (66) 



By collecting all spatial components we get the contribution ( I49ll corresponding to a 
particular cycle 



d"^giro((Pq),q;t) 



-j n 2 d^'qi 
n 



m 



271 hit 



N^'Vd, (67) 



where we reintroduced the notations (jUD and n = N^^ for the length of the particular 
cycle under investigation. Note that this general form also includes the case of a one- 
cycle A'';^ = 1. 
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By considering all traces corresponding to the cycles of one particular permutation 
one gets 



d^^giro((Pq),q;t) 



m 



2Tihit 



(68) 



U)=l 



(16811 as an expression associated with a permutation P only depends on the partition of 
into cycle lengths (as does the measure ii{M. p) fj44|) ). By collecting all permutations with 
the same partition in the sum over Sn-, the trace of the symmetry projected propagator 
can be written as 



N N 

/ d^^g iro,±(q, q; t) =— $^(±1)^^' '^^-^iv. c(iVi, 

' 1=1 Ni,...M=l 



Ni) 



Ni,....Ni=l 
Ni<---<Ni 

ID I 

/ m \ — 



(69) 



UJ=1 



As before, c{Ni, . . . , Ni) denotes the number of permutations with a cycle decomposition 
of lengths Ni, . . . ,Ni fH^ . Using the bilateral Laplace transformation rule 



-1 



rr(^)l 










2m J,^ 



z/ > 0, 



(70) 



with the Heaviside step function 6{x), the final result for the smooth part of the 
symmetry projected DOS f|T9|) for a system of identical non-interacting bosons (+) or 
fermions (— ) in a D-dimensional billiard of volume Vd without confinement corrections 
reads 



1 ^ 

1=1 



,N-l 



N 



iVi,..,Ar,=i 

Ni<-<Ni 



' 1 

LtJ = l 



(71) 



X 



m 



2 V ] 



D 



2nhy r(f) 



E--'e{E) 



In general (17T|) is a sum of powers of E with coefficients that are, besides their 
dependence on the billiard volume, expressed as sums over partitions N = Ni + ■ ■ ■ + Ni 
only depending on A^, / and D and therefore universal for all A^-particle billiard systems 
with dimension D. A more compact form of ( 17Ti) can be obtained by appropriately 
scaling the density and energy and rewriting the sum over partitions as ordered tuples. 
Using (H6i) and defining the scaling- density 
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which is roughly the density at the first single-particle level, and the universal coefficients 

^ ' 1 ^4-1 



leads to 



Ni,...,Ni=l uj=l 



SJE) = f (±ipc, <^^9(i5) . (74) 

An extension of fl71|) is given by the inclusion of physical boundary effects. Its 
derivation under the assumption of locally flat boundaries can be performed directly 
by incorporating the additional geometrical features to the propagation in cluster 
zones. But as this way to proceed gets extensive and seems a bit long-winded, an 
alternative, indirect but equivalent derivation utilizing a convolution formula fl82|) by 



Weidenmiiller |2T] is chosen for this publication and can be found in Appendix C The 
resulting expression reads 



UE) = ,0 E(±i)"" E TT7T %lrr ' 

where A = l\D/2 + ls{D — l)/2. It depends on the universal coefficients 



Ni,...,Ni=l iu=l w=/v+l 



and the dimensionless geometrical parameter 



l = ±^^, (77) 

D 



4:?/V, 



which represents the ratio of the surface Sd-i to the volume Vd of the billiard. The 
plus(minus) sign in f l77|) refers to van Neumann(Dirichlet) conditions at the physical 
boundary. Equations fl74|) and fl75|) can be regarded as the main results of this 
publication. ly + ls = l represents the splitting of all / clusters into such that contribute 
via free propagation in the interior (/y) and such that contribute by reflection along the 
boundary (Is) of the billiard. In the case of D = 1 the summands corresponding to 
Zv = have to be replaced according to the rule 



r(M£^) 



eiE) ^ 6igoE) (7^ 



and in (177|) the surface 5*0 has to be taken as the number of bordering end-points, which 
would be two for a finite line and zero for the one-sphere-topology. The unconfined 
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expression (171|) can easily be reobtained as the special case 7 = by recognizing that 



The highest power of the energy in ( 17^11751) has the exponent {ND)/2 — 1, which 
reminds of the Thomas- Fermi approximation of the effectively iV-Z)-dimensional billiard. 
This is not surprising, since / = N refers to a partition into unities = 1 + ■ ■ ■ + 1 
associated solely with the identity permutation P = id^^ = ( 1 ) ( 2 ) ■ ■ ■ ( ). In the 
geometrical picture, this corresponds to the propagation of individual particles. None 
of them are clustered and in ( 1751) none of them are reflected on the boundary, since 
/s = for the highest power. The combinatorial factor (146|) and the coefficients (!73ll 
and ( |76l) compute to c(l, ...,!) = Cn = Cn,n = 1 and the corresponding term in (|7T1) . 
([71]) and fl75l) is the volume Weyl term of the fundamental domain in ■ D-dimensional 
space 



The next section will show the importance of all corrections in fITT]) . fITl]) and (1751) 
beyond the volume term (!79|) . 

8. The Geometrical Emergence of Ground State Energies in 
Fermionic Systems 

As often in the context of semiclassics the case of a two-dimensional billiard is of special 
interest. On the one hand, this is because of possible technical applications. One can 
think of confined two-dimensional electron gases in semiconductor heterostructures or 
two-dimensional superconducting structures with bosonic description due to Cooper 
pairing for example. On the other hand, the existence of equally distributed energies in 
a 2D single-particle billiard without boundary corrections is a valuable special feature. 
This is not only because of the exceptionally simple form that the density of states takes 
in these systems. A constant single-particle smooth part also opens the possibility to 
make connections to number theory. Namely approximations for average distributions 
of partitions of integers can be related. 

For positive arguments E the unconfined DOS ( 17^ in Z) = 2 is a polynomial 
of degree A^ — 1 in the energy with coefficients that are just rational numbers. The 
coefficients can be summed up exactly for explicit values of A^ and / albeit with 
computation time increasing very strongly with A^ when using the form at hand fl7Tl) 
or ( 17^ . Note that the form of the coefficients used in ( 17Ti) in terms of ordered partitions 
corresponds to less computation time while the form (17^ seems to be a better starting 
point for simplifications or analytical calculations. 



Ci^i^=i — Ci. 




(79) 
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Figure 5: Symmetry-projected DOS without boundary corrections m. D = 2 for bosons 
(green, dotted) and fermions (blue, dashed) in comparison to the naive volume term 
(red, solid) for several numbers of particles. The vertical dashed black line shows the 
particular expected many-body ground state energy for fermions (|80l) . Densities q 
are measured in units of the constant single-particle density Qq. Energies E are measured 
in units of its reciprocal q^^ 



Figure [5a] shows the case of two particles. The bosonic and fermionic cases are shown 
in comparison to the naive volume term 0791) . Already here, the symmetry corrections 
give qualitatively the right picture. With respect to the naive term, the fermionic 
density is shifted to higher energies, which is according to the expectation of the many- 
body ground state energy below which effectively no level should appear. Whereas the 
bosonic density is shifted to lower energies, which accords to the full counting of many- 
body levels corresponding to shared single-particle energies in contrast to the naive 
term, where these are counted with a factor of even if they can not be permuted 

in N\ ways due to identity of some of the single-particle energies. 

Figures [5b] - [5fl show the cases of two to twenty particles. In the fermionic case 
the lower powers in E in the polynomial produce oscillations around the axis ^ = 0. 
With increasing particle number these oscillations get smaller in amplitude and larger 
in number, approximating a zero-valued DOS. The DOS is effectively shifted to higher 
energies and an energy gap opens that coincides with the expected fermionic ground 
state energy i?^g calculated by counting single-particle levels by virtue of the smooth 
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Figure 6: Symmetry-projected 
DOS m. D = 2 for bosons (green, 
dotted) and fermions (blue, dashed) 
in comparison to the naive volume 
term (red, solid) for = 20 
particles. Negative values of the 
fermionic density are plotted as 
logarithm of its absolute value. The 
expected many-body ground state 
energy for fermions flHOl) is 

represented by the dashed vertical 
line. 



single-particle DOS Qsp{E). It is defined by 



(f) .. 

GS •" 



dE' -Q,^{E')E' . 



(80) 



where the Fermi energy E-p is defined through 



dE'0,p{E'), 



^1) 



but instead of explicitly filling up single-particle energy levels by hand, this time the 
ground state energy occurs as a consequence out of exchange symmetry incorporated as 
a modification of the propagator. The corrections from cluster zone propagations are 
sufficient to automatically reproduce the expected ground state energy. When increasing 
N, the symmetry projected DOS at = E^^ is observed to keep moderate values 



"'(f)^ ~ 0{1) while the naive density at this energy grows exponentially with A^. In 
contrast to the fermionic density, the bosonic density does not have these oscillations. 
There, the polynomial in E has only positive coefficients and the density is effectively 
shifted to lower energies as expected intuitively. 

For higher particle numbers a single logarithmic plot suggests itself in order to show 
at the same time the oscillations that are becoming very small and the strong growth 
behaviour around and above the ground state energy. Figure [6] shows the smooth part of 
the density in the case of iV = 20 particles. Again the fermionic energy gap accurately 
reproduces the ground state energy, indicated by crossing the axis of abscissa. Also in 
the bosonic case, the corresponding density g+{E) apparently keeps moderate values at 
the expected bosonic ground state energy Eg] computed in analogue manner to 4^ 

The small values of Q-{E) for E < i^^g result from large cancellations regarding 
the different terms in the sum fl74l) with different values of /. Therefore the behaviour 
of the DOS in this regime is very sensitive to numerical errors and all the corrections 
are needed to reproduce it correctly. In order to illustrate this fact figure [7a] shows the 
deviations one obtains when leaving out contributions. 
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Another feature of the oscillations around zero np to E = Eq^ is their rigidity with 
respect to integration in the sense that not only the DOS itself oscillates around zero 
but do also several integrals of it. Some integrals are shown in figure [Tbl 




0.001 - 



(a) Absolute value of the many-body DOS for N — 
13 fermions. The blue solid curve shows the full 
expression (|7ip . The black curves correspond to 
leaving out the contributions I = 1 (dashed), I = 2 
(dash-dotted) and / = 3 (dotted). 




(b) Absolute value of multiple integrals 
g(-"\E) = de„ • • • dei e_(ei) of the 
many-body DOS for N = 17 fermions. 
Oscillations around zero for energies below 
^GS (vertical dashed line) remain. 



Figure 7: Sensitivity (figure [7a|) respectively rigidity (figure [7b|) of the oscillations of 



The effect of boundary corrections is shown in figure [HI Displayed is the = 12 
fermionic unconfined level counting function Af{E) = dE' (£") in 2D and in 
addition to that the corresponding confined case (1751) for Dirichlet boundary conditions 
with a geometrical perimeter-to- area-ratio of 7 = — -\/7r/2, which is the smallest possible 
parameter for singly connected billiards since it is the one of a circular billiard. The 

curve seems to be shifted further to higher energies, enlarging the energy gap. Again the 

( f) 

expected ground state energy E^^, this time calculated using the single-particle Weyl 
expansion with perimeter correction, is reproduced very well. Already for this minimal 
value of 7 the deviations from the unconfined case apparently are rather strong. Thus 
naturally the question arises whether the assumption of locally flat boundaries gives a 
sufficient description of an actual billiard or additional corrections also lead to rather 
strong deviations, making the treatment of curvature inevitable. Therefore the exact 
quantum mechanically solvable levels of a circular billiard have been arranged to non- 
interacting fermionic many-body levels shown by the green staircase function. The 
deviation seems again to be effectively a slight shift, smaller than the deviation of the 
confined from the unconfined case. Not displayed is E^^ including curvature, which in 
this case coincides with the exact circle levels, in the sense that it lies half way between 
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N = 12 




E[Qo] 

Figure 8: Level counting function for 12 fermions. The blue dashed curve shows the 

smooth part without confinement corrections. The black dotted curve includes boundary 

corrections with a geometrical parameter of 7 = — -\/7r/2, which corresponds to a circular 

billiard. The green staircase (left) shows the exact non-interacting many-body levels of 

the circular billiard. The red staircase (right) shows the exact levels of a cylindrical 

billiard with same geometrical parameter. The dash-dotted black curve in the inset 

( f) 

shows the smooth part shifted to lower energies corresponding to the shift of Eq^ due 
to curvature. 



the first two of them. The additional comparison with exactly computed levels of a 
billiard with the shape of a cylinder barrel with same 7-ratio as the circle serving as 
an example of a 2D system without curvature shows good agreement with the smooth 
part. Note that in this geometry the minimal value of 7 could be underrun, which is not 
contradictory since it is not singly connected. The absence of curvature can be seen in 
the single-particle expansion given by Balian and Bloch [26], where the corresponding 
correction is proportional to the Euler characteristic x of billiard, which happens to 
be one for a disk and zero for a cylinder barrel. 

An analysis of the Eq^ in 2D for the unconfined and confined cases with and 
without curvature indicate the relative importance of the corresponding contributions in 
the smooth DOS. The smooth ground state energy involving a Dirichlet type perimeter 
correction without curvature (x = 0) reads 

2 ]. II 

4's(^, 7, 0) = 4's(^, 0, 0) (v^TT^ + v^) {VTT^ - - v^) , a = . 

The inclusion of the curvature correction x^{E)/Q [26] in the single-particle DOS yields 
the ground state energy ii^Qg(iV, 7, x) = Eq},{N — x/6,7,0). Comparing the correction 
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Figure 9: Level counting function of a ring billiard with = 6, 9, 12, 16 fermions. 
The blue dashed curve shows the (absolute value of the) smooth part with confinement 
corrections ( !75|) using a geometrical parameter of 7 = — -\/37r/2. The red staircase shows 
exact non-interacting many-body levels of a ring billiard with the same geometrical 
parameter which corresponds to a ratio of radii of 2. The black dotted curve displays 
the Bethe approximation ( 1103P using the exact ground state energy as input. The 
horizontal black line corresponds to A/" = 1/2. 



due to the perimeter = 0) with the further correction due to curvature 

Eg{N,^,0)-Eg{N,0,0) 8|7|ViV 16iV ^ ^ 

indicates that curvature contributions in general get strongly suppressed for large 
particle numbers. Based on this argument, the smooth part of the many-body DOS 
including curvature in D = 2 might be approximated by simply shifting the many- 
body DOS (including perimeter corrections) by AE = EQg{N, 7, x) ~Eqs{N, 7, 0). The 
corresponding function is plotted in the inset of figure [81 It has to be emphasized that the 
ground state energies reproduced by the smooth densities are defined by smooth counting 
and do not necessarily have to accurately coincide with the exact many-body ground 
state energies of actual non-interacting quantum billiards. The exact fermionic ground 
state levels usually are subject to fluctuations around Eq^ with respect to the particle 
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number A^, which can be treated with semiclassical methods including corrections from 
periodic orbits [30]. For rather regular systems possibly featuring degeneracies, these 
fluctuations in general can become strong and even the average (£'Qs(A^))iv over some 
window in N can deviate from the smooth prediction Eq}.{N) [30]. A related feature 
of exact non-interacting many-body spectra are fluctuations of the many-body DOS at 
low excitation respectively temperature which can also be reproduced as semiclassical 
shell corrections to the Bethe estimate fll03p as shown by Leboeuf et al [31]. Therefore, 
the smooth part of the symmetry projected DOS ( 17411751) presented in this work should 
not be regarded as an estimate for the exact ground state energy of a given system 
or the behaviour of the actual many-body DOS at energies close to it. Nevertheless, 
as fluctuations tend to wash out for higher excitations, the correct behaviour is very 
well reproduced for energies above a critical excitation (which depends on the specific 
geometry of the system and the particle number N). Note that the theoretical prediction 
of this more universal behaviour does not need any further system specific input than the 
volume and surface (and a possible non-zero valued Euler characteristic) of the billiard. 

Figure I9] shows both, deviations of the smooth prediction near the ground state 
as manifestation of such fiuctuations and the increasingly good agreement for higher 
energies using the example of a ring billiard with several numbers of fermions. Like 
the cylindrical billiard, this system has an Euler characteristic of zero and therefore 
there are no signatures of curvature. Compared are the smooth level counting function 
M{E) in the confined case and the exact many-body counting function of a ring billiard 
with a ratio of i?/r = 2 of outer and inner radius, which corresponds to a geometrical 
parameter of 7 = —\f^/2. 

9. Comparison with Known Results 

9.1. Smoothing in Weidenmiiller Convolution Formula 

As Weidenmiiller pointed out [2Ij, the symmetry projected part of the exact DOS of a 
non-interacting system can be constructed as a sum of convolutions of the single-particle 
DOS with modified energies. The DOS for a non-interacting bosonic (fermionic) system 
with exact single-particle DOS &p(-E) = Yli^i-^ " ^i) reads 



^ N N I 

^±(^)=mB±1)''~' E ^^,E^. c(Ar„...,Ar,)(n^) 




X 




fdE,...dE,6[E-J2^-) ri^^pd^) 

,.,—1 ,..—1 ^ 



The derivation of fl82|) is on the one hand based on the fact that the non-interacting 
many-body propagator decomposes into a product of single-particle propagators. On 
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the other hand, this form requires the convolution property 



/ 



dVir(^p)(q",q';t2)ir(^P)(q',q;ti) = TT^^p) (q", q; ti + ts) 



(83) 



of the single-particle propagator, which is strictly fulfilled for the exact quantum 
propagator associated with the exact DOS. Inserting the single-particle spectrum 
directly into (l82l) . one can, in accordance to [21] show that (182!) is a formal way to express 
the construction of many-body levels out of filling up single-particle levels under the 
condition that these should be counted correctly. So that no two filled fermionic single 
particle levels are the same and bosonic states where any of the filled single-particle 
levels are the same are not undercounted by the factor but instead fully addressed 
by counting each of the corresponding possible distributions of single-particle energies 
among the N particles with the factor of (Yloj ^lu^-) I Here, the N^^ denote the numbers 
of particles sharing a single-particle level in the questioned many-body state. 

The question of whether and, if yes, to what extent formula ( 182|) can be utilized 
in order to obtain or compare to the smooth part of the symmetry projected DOS can 
be answered in a twofold way. Naturally this question reduces to the question of how 
reasonable it is to simply use the smooth part of the single-particle DOS in fl82|) instead 
of the exact one. 

First, one could regard the smooth part of a spectrum as a convolution with some 
smearing function. The question is then, if the smoothing of the single-particle DOS 
in (1821) yields the same result as if one would smooth the many-body DOS directly by 
a similar smoothing procedure. It turns out that this can be answered positively if 
the function used for smoothing obeys two conditions. Let 6a{x) denote the smoothing 
function with a continuous parameter a describing its sharpness and J^^dx 6a{x) = 1. 
If this function fulfils 



then the replacement of the single-particle DOS Qsp by its convolution with 6a yields the 
same as convolving the many-body DOS g±. with S^a- If the function 6a perceived as a 
distribution had a finite mean value and variance, then implies that a is proportional 
to the standard deviation. But then this would contradict the second condition (185|) as 
the variance in general is additive under convolution instead of the standard deviation. 
This shows that a proper smoothing function in this context should not have a finite 
standard deviation. Nevertheless, there is a distribution without standard deviation 
fulfilling the requirements at hand, namely the Cauchy distribution 







(86) 



The second way to address the question is to replace the exact single-particle DOS 
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in ([H2D by the smooth part ^sp defined by the Weyl expansion up to a specified order. 
In the cases of restricting to the volume term and considering also the first boundary 
correction in the corresponding integrals over single-particle energies are carried out 
successively yielding a solvable recursion relation. [Appendix A| and [Appendix C| show 



these calculations with the outcome that this procedure yields the same result for Q±{E) 
as the investigation of the propagation in cluster zones fITil [75]) . 

The reason why both procedures give the same result lies in the requirement f[83|) . 
For the exact quantum propagator the fulfilment of it is obvious. In the unconfined 
case the exact propagator is replaced by the free propagator, which also obeys the 
convolution property, but only if the intermediate coordinates q' are integrated over full 
space and not only the interior Vl of the billiard. Recalling the derivation of ( 17^ . 
this corresponds exactly to the extension of the integration limits for all intermediate 
coordinates g2, •••,<?« in ( |60i) allowing the evaluation of the integral as multidimensional 
Gaussian. Thus the assumption of infinite distance to the boundary in the analysis 
of cluster propagations in combination with the usage of the free propagator actually 
preserves the convolution property f[83|) . 

In the confined case the free propagator is modified by a refiection term 

Kt\^. q; t) ^ ^^^(q', q; t) ± <"^(i?q', q; t) , (87) 

where R denotes the refiection with respect to the boundary regarded as locally fiat and 
± refers to Neumann respectively Dirichlet conditions. This replacement yields the Weyl 



expansion including the surface correction. The analysis in [Appendix B[ shows that the 
so constructed propagator in combination with the assumption of local fiatness of the 
boundary also exhibits the convolution property ( l83i) . This shows the equivalence of the 
computation Q±{E) in the confined case (!75|) via the convolution of single-particle Weyl 
expansions up to the first boundary correction and the computation by investigation 
of the propagation in cluster zones including refiections on the boundary assumed as 
locally fiat. 

This equivalence gives rise to the question whether it is possible to relate corrections 
from propagation in cluster zones to the correction of delta peaks that is inherent to the 
exact convolution formula ( |82l) of Weidenmiiller. And indeed one can relate each cluster 
zone correction to the correction of delta peaks for total energies. The corrected total 
energy is the energy that is a partition of single-particle energies just the way the cluster 
zone corresponds to a partition of all particles into clusters. The cluster zone correction 
is associated to the correction of the total energy that is built of single-particle energies, 
where all particles in a cluster share the same energy. But of course it is doing it in 
a smooth manner, meaning it produces a smooth correction to the DOS corresponding 
to a smooth version of the density constructed out of such many-body energies with 
correct pref actor. 
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9.2. Connection to Bethe's Estimate in Fermionic Systems 

This section is dedicated to the connection between the expressions presented in this 
work and the celebrated Bethe approximation to the many-body DOS. As the latter 
results as a saddle point approximation of an integral representation of the DOS obtained 
by general quantum statistical methods for non-interacting fermionic systems, the link 
between both is the equivalence of the smooth part of the DOS (17411 75 p and the standard 
statistical formulation in terms of the single-particle DOS. This connection will be drawn 
by the usage of the technique of generating functions as the easiest to handle statistical 
object associated with the DOS is the grand canonical partition function Z'^'^^z^P) 
which itself is the generating function of the canonical partition functions for N particles 



ZC(iV,/3) 



1 d 



N 



^^^(^,/3) 



A^! d^^ 



2=0 



(88) 



The essential step in order to find the generating function for the DOS ([75 



G±{z,E) = Y,Q±{N,E), 



,N 



^Q±{N, E) 



N=0 
1 d 



TV 



A^! d^^ 



G{z,E) 



2=0 



(89) 



(90) 



is to represent the difficult to handle combinatorial numbers Ci, Ci^i^ ( ITTl) . ( 1751) in terms 
of their generating functions. In order to ease notation, the abbreviations fi = ^ + 1 



and u = + 1 are used. 



and similarly 
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(TV) 



N 



1 d 



^n=0 



TV 
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(91) 



(TV) 
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W\Az^ 



E- E- 



n=0 



n=0 



TV 



1 d 

W\dz^ 



[{U,{z)t {U^t 



2=0 



(92) 



2 = 



These relations can be seen by factoring out the powers of the sums in parentheses 
and recognizing that the A^-th power of z in the product must be built of / factors 
i = 1, . . . , / with Ui = N. Thus each such combination corresponds to a particular 
partition of A^ into / parts ni + ■ ■ ■ + n/. The correct coefficients are then created 
as products of the parts n'^^'^ appearing as denominators in the sums over n. The 
generating functions can also be built constructively which will be shown exemplarily 
for the unconfined case. The confined case can be treated similarly. For this calculation 
the form of C^^^ expressed as sum over distinct partitions instead of ordered tuples will 
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be used. Note that C^'^^ should be understood to be zero for all / > A^. The starting 
point is to calculate constructively the generating function of C\^'^ /l\ which is 



9i{z) 



N=0 N=0 Ni<...<Ni=l V=l / ^uj=l 



E n 



(93) 



where the sum over all possible distinct partitions is expressed as sum over all 
multiplicities m„ with which the part sizes n = 1,2, .. . appear. Given a fixed value 
of / corresponds to the restriction = I in the sum. This restriction is dropped by 

constructing a further generating function of fl93p . this time with respect to /. 



E 



9i[z)y 



=1 



n 



n=l 



y— 



1 / z"-y 



m=0 



ml \ 



n 



n=l 



exp 



z y 



exp [yUf,{z)] . 



(94) 



Making use of this secondary generating function shows that 

. . Id' 



n dy 



7 exp [yU^iz)] 



-0 /! 



(95) 



Note that because of Li^(2;) = C(z) as 2 — )■ 0, indeed the Ci vanish for / > A^. 



From here, the confined case will explicitly be treated, including the unconfined case as 
7 = 0. Using (I92D in (CSl) gives 



lv,ls=0 



r(A) 



(96) 



with the abbreviated exponent A = l\D/2 + ls{D — l)/2 and /i = D/2 + 1, 
V = (D — l)/2 + l as before. The upper limit of the sum over / has been raised from 
to infinity which is no problem since the A^-th derivative at z = is not affected by 
terms with / > A^. Applying Laplace transformation switches from energy domain to 
the domain of inverse temperature and allows for exact summation of /y and /s- 



G±{z,E) = C',^ 



E 



[±U,i±z)y-[±U,{±z)y^ ^ / go 



7 



/3 



exp 



f±mX(±.)±7f'° 



V/3 



UJ±z] 



(97) 



In the case of D = 1, the replacement rule ( 1781) is already accounted for in (1971) . Note that 
the bilateral definition of Laplace transform is used, which simply yields the Heaviside 
step function 6{E) after taking the inverse Laplace transform. 
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The result (l^7|) is equivalent to the grand canonical potential one obtains by 
standard statistical methods. In general, the grand canonical potential of a non- 
interacting quantum system of identical particles can be expressed in terms of the 
single-particle DOS as 

lnZ^C(z,/3) = T J dEg,^{E)HlTze-^^), (98) 

where the upper sign corresponds to bosonic symmetry {z has to be small enough in 
modulus to keep the argument of the logarithm inside the unit disk around 1) and the 
lower one corresponds to fermionic statistics. If the single-particle DOS is taken as its 
smooth part built as a sum of smooth functions of the form 

Q^iE) = ^o7.^^Jp^(i?) , (99) 
each of the summands yields an additive contribution 

lnZ«J(.,/3) = ±7.^.^'f;^^-^ r dE E--'e-''^^ 



oo / , xfc 



' k=l 

/3 



Li.,+i(±z), (100) 



which shows the equivalence to (jnZD due to Z*-^^ = C[G] when taking the appropriate 
single-particle DOS with 71 = 1,1/1 = D/2 in the first term and 72 = 7, z/2 = (-D — l)/2 
in the second one. The derivation of (198 p involves the replacement 

5^/(E,)-^ f dE&p(E)/(E) (101) 

i 

of the sum over single-particle eigenenergies, which are assumed to be discrete in 
the first place. Simply taking the smooth part of the single-particle DOS instead is 
therefore closely related to taking the smooth part in the convolution formula (152]) 
by Weidenmiiller, which we showed to be equivalent to the actual smooth part of the 
many-body density obtained with the geometrical approach. 

For practical reasons it should be mentioned that in the unconfined case 7 = in 
D = 2, the inverse Laplace transform in ( |9711 can be done explicitly to get 



,(D=2,7=0)/ _ ^ / , Li2(±^ 




G'^ '' >{z,E) = g,^l±^^ l,(2^±U,{±z)goE) , (102) 

which could also be seen directly from f PB]) after recognizing the power series of the 
modified Bessel function Ii(x) with the two factorials lyl and r(A) = (/y ~ 1)' in the 
denominator (/g = 0). 
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The Bethe approximation [13] of the non-interacting fermionic many-body DOS is based 
on the general statistical relation ( l98l) and reads 



where Q = E — Eqs > is the excitation energy above the many-body ground state 
energy (15U]1 and g = gsp{Ep) is the single-particle DOS at the Fermi energy ( 15T]) . It 
results from a low temperature expansion and a simultaneous complex saddle point 
approximation in two integrals. The first one being the Bromwich integral achieving 
the inverse Laplace transform of and the second one being a closed complex 

contour integral representing its A^-th derivative with respect to z. The low temperature 
expansion limits the validity of f ll03p to low excitation energies Q < Ep (range increases 
with A^) whereas the defect produced by the saddle point approximations becomes 
noticeable when approaching the ground state energy Q ^ g^^ and culminates in the 
divergence a.t Q = 0. 

Despite the fact that (I103P is valid in general and is not depending on the specific form 
of the single-particle density, it shall be sketched here how to obtain it explicitly for the 
unconfined case using expression fl971) respectively fllOOp . since the explicit expressions 
can be of use investigating possible extensions to f ll03p that are not general. Expressing 
the inverse Laplace transform and derivation with respect to z as integrals, the smooth 
part of the fermionic many-body DOS for particles in a D-dimensional billiard 
(z/ = D/2) without confinement corrections reads 

g_{N, E) = (—^ [dp [da exp [^E - aN - /3-'^Li^+i(-e")] , (104) 
\27nJ Jy^ Jr^ 

where e" := ^ and energies are measured in units of ^o- The contours are F/j : 
(e — ioo, e + ioo) with e — )■ 0"*" because of the essential singularity at /3 = and 
Fq : (a — in, a + iir) with real a < 0, corresponding to a closed counterclockwise circular 
contour for z inside the unit circle, chosen because of the branch cut z G {—oo, —1] of 
Li;^+i(— 2;). Applying saddle point approximation in both integrals will yield a saddle 
point that has a large positive real value in a in the requested regime. Therefore 
the contour for 2; = e" has to be thought of to be deformed outside the unit circle 
enclosing part of the branch cut. Nevertheless, only the behaviour in the vicinity of 
the saddle point shall be regarded, dropping the exact form of integration contours and 
neglecting the integration along the branch cut, allowing for an asymptotic expansion 
of the polylogarithm for large arguments \z\ ^ 1 or 3?(a) ^ 1. The way to proceed 
is the following. First a is considered to have a large real part. From the consequent 
approximation the saddle point will be found to fulfil this statement, justifying the 
assumption afterwards in a self-consistent manner. The asymptotic expansion of Li^+i 




(103) 
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leaves the exponent in (11041) - which has the statistical interpretation of entropy - as 



The saddle point equations are 

at the saddle point values of [3 and yU = statistically interpreted as inverse 

temperature and chemical potential. The right hand sides of (11061) and (11071) match the 
statistical definitions of mean particle number (A^) and mean energy {E) in the grand 
canonical formalism. Thus the saddle point values of /3 and /i fix the ensemble averages 
of energy and number of particles to the given values (A^) = and {E) = E. For 
further computation, it is convenient to express the equations in terms of the following 
quantities. 



rw r(. !)■ 

with the level counting function Af{fi) = dE g{E) and the total many-body energy 
S{fi) = J^^dE Eg{E) up to = /i. ^sp has been abbreviated to g in order to ease 
notation. The saddle point equations then read 

^2 



TT 



N = Af{fi) + —^g\fi), (109) 



E = S{f,) + — {f,g{fi)y, (110) 

matching the Sommerfeld expansions up to the first term. Since fi ~ Ep at the saddle 
point in the requested regime, one can expand 

^ AfiEp) +g{fi) (/i - Ef) , (111) 

=N 

=> S{fi)^S{EF)+fig{fi){fi-EF) ^ EGs-^/if?'(/i), (112) 



leading to 



6/32 

The entropy at the saddle point 



^(/i). (113) 



S = (l + ^^(3E ~ aN (114) 
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13 ' 'r(z/ + 2) ' '6/32 r(z/) r(z/ + i) 6/32r(i^-i)' 



And by exploiting flllSp one finds 5 = y -|- gQ. What is left is the computation of the 
determinant of second derivatives of S with respect to /3 and a 



det S"\ 



dmAE)) 
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dmAE)) 


d{a,f3) 




difi,f3) 



(116) 



After neglecting terms that are sub-dominant in the low temperature limit like 
i < g^'^Kf^), one finds 

I det ^"1 ^ ^[^(/x)]^ = - • 

Collecting everything and accounting for a minus sign due to the directions of the 
integration paths in the saddle point finally gives the Bethe approximation 



g_{N,E) ^ ^=-^--—exp { xl^g ■ (E - Egs) ] , (118) 



/48(E-Egs) 

with the ground state energy Egs = Jy[T{iy + 1)NY+-/T{iy + 2) which matches the 
definition of smoothly counted ground state energy E^g (IHOj) . What is left is to show 
the consistency of « ^ 1 at the saddle point. From 01131) one obtains 

^ _ 7r/ig(/i) ^ nuMjfi) 

Together with = N — -^^{^ — l)A/'(/i) from (11 111) , a computes to 

a^-==--^ L. (120) 

ybgQ D a 

Thus in the regime of low excitation energies and high numbers of particles, the 
assumption is self-consistently justified with 

a^^ + 0(^). (121) 



Figure [10] compares the Bethe approximation with the unconfined densities in two 
dimensions for different numbers of particles. The approximation gets better for 
increasing and decreasing Q (down to Q ~ Qq^)- The displayed energy window 
scales with A^^. Thus the similar looking plots for A^ > 10 show that the range of 
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Figure 10: Symmetry-projected DOS without boundary corrections in D = 2 for 
= 3, 5, 10 and 28 fermions (blue, solid) in comparison with the Bethe approximation 
(black, dotted) and the finite-size corrected improved form fll22p (brown, dashed). 
Densities are measured in units of dZS]). The grey line indicates g{E) = qq. 



excitation energy in which QsiQ) gives reasonable results roughly scales also with A^^ 
as indicated by fll2ip . Note that the Bethe approximation needs the information about 
the particular many-body ground state energy as additional input (which has here been 
taken as Eq^) whereas the smooth part of the unconfined symmetry-projected density 
does not but instead automatically provides it. In addition to QbIQ), also an improved 
form including finite- iV corrections. 



gF{N,Q) = QsiQ) exp 



(122) 



is plotted. Equation fll22p can be found in number theoretical context [32j as 
approximation to the smooth part PN{n) of the number of possible partitions of an 
integer n where all parts are restricted to a maximum size N. This function is equivalent 
to the non-interacting many-body DOS for an equidistant single-particle spectrum which 
makes the connection between number theory and many-body physics in this context. 
There is also a generalization [22] of (I122p obtained by saddle point approximation in 
the statistical approach similar to the derivation of Qb{Q) but including a correction 
~ e~" to the asymptotic expansion of the polylogarithm fllOSp in fll04p for small ^{a). 
The resulting expressions are generalized for single-particle DOS of arbitrary power 
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10. Conclusion and Outlook 

In this paper we have made explicit the essential geometrical content behind the smooth 
part of the density of states in systems of non-interacting particles. We have shown that 
both the functional form of the level density and the ground state energy are consequence 
of a large cancellation effect between contributions polynomial in the energy related with 
manifolds in configuration space which are invariant under the group of permutations. 
Our results cover the whole regime of energies and are valid for any number of particles. 

The main message of our work is the rigorous construction of the connection 
between polynomial contributions to the density of states and the measures and 
dimensions of manifolds in the classical phase space which are invariant under particular 
permutations of the particles. Moreover, by means of algebraic manipulations we have 
transformed our expression in a way that explicitly shows its equivalence with several 
known but previously unconnected results. In particular, by re-expressing the sum over 
cluster zones in terms of generating functions, we have established its equivalence with 
the thermodynamic formalism leading to the celebrated Bethe estimate for the level 
density in fermionic systems. 

The geometry of invariant manifolds presented here is a purely kinematic concept 
independent on any particular dynamics, and although we have used it to fully construct 
the level density for non-interacting systems fl74|75p . we expect that it plays also a 
fundamental role when interparticle interactions are present. In order to apply our 
formalism in this more general scenario, the key step is to relax the condition expressed 
by (l47ll which makes explicit the assumption of non-interacting particles. In the case 
of interacting particles, however, the key concept remains the same: the smooth part 
of the density of states does not require the solution of the problem including both 
interactions and confinement simultaneously. Symmetry effects on the Weyl expansion 
for interacting systems are encoded in unconfined (but interacting) propagation around 
cluster zones. 

In order to keep the problem tractable, well controlled approximations for the 
unconfined (but still interacting) many-body propagator are necessary, and the first 
task to extend our methods into interacting systems is to check the physical picture one 
gets when a certain approximation for the unconfined propagator is used around the 
cluster zones. One possible scenario where this program can be carried on is when we 
treat the interactions in eikonal approximation, freezing the interactions between sets 
of particles in different cluster zones. This and other possible approximations will be 
reported in a future publication. 
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Appendix A. Calculation via Convolution Formula - unconfined case 

The single-particle Weyl volume term reads 



Using ( 1821) yields the sum 
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(A.4) 



UJ = 1 



X{E) can be calculated by successively performing the integrals of single-particle energies 
E^ and solving an emerging recursion relation. It is convenient to write r = D/2 — 1 
and define a new variable for the energy to be distributed among the first n particles 
for every n 



n = / - 1. 



(A5) 



The first step in the integral of I{E) is then 



dEiE[9{ei - El] 
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(A.6) 



The next integral over E2 is of the form 



\x £ {c - xY9{c) = c-+^+i£(i±l)£(i±^^(c) 



V(2 + r + s) 



(A.7) 



with c = 62 = e^ — E^. Therefore, also the third integral and all others are of this form. 
Let An be the total prefactor and Sn the exponent s in ( ]A.7I) appearing in the n-th 
integral step. Then one gets the recurrence relation 



A ,-A ni+r)ni + sn) 

Sn+l = r + Sn + l, Si = 0. 



(A.8) 
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The convolution (1A.3P is then expressed as 

^(^) = (ri W') ^ ' E^'-^-'e{E) . (A.9) 



U}=1 



The solution of the recurrence reinserting r = D/2 — 1 is 

Using the the original expression for a one gets the final result 

^conv,±(i?)=^5^(±ir-' 5^ 5^,j:;v. c(Ari,...,ArO 

i=l Ari,...,Af,=l 

7Vi<...<JV, (A.ll) 

ID 1 



X 



tJ=l \ ^ / 



which is identical to the result obtained by investigating the propagation in cluster 
zones (TfTj) . The usage of the volume term in the single-particle Weyl expansion is 
consistent with neglecting the single-particle billiard boundaries in the calculation of 
short path propagations. 



Appendix B. Convolution Property of the Confined Propagator 

Let the confined propagator for a single particle in first correction be denoted by 

K{ci2, qi; t) = Ko{q2, qi; t) ± Ko{Rq2, qi; t) , (B.l) 

where ± refers to van Neumann or Dirichlet-boundary conditions respectively. Rq2 
denotes thereby the reflection of the point q2 with respect to the boundary, locally 
regarded as being fiat. In the spacial integral that has finally to be performed to obtain 
the DOS, for points lying far inside the billiard compared to the wavelength in question 
the second term won't contribute. This justifies regarding R as the unambiguous 
reflection with respect to the tangent plane in the nearest boundary point. Thus, 
whenever a point far inside the billiard is formally reflected in the expressions below, 
one should keep in mind that then the expression won't contribute and that therefore 
any ambiguity does no harm. The assumption of flatness also implies 

iro(q',q;t) = i^o(i?q,i?q;t), (B.2) 

where q', q are arbitrary positions not restricted to the interior of the billiard. The 
convolution of two propagators where all of the three involved positions (initial, final 
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and intermediate position) are lying in the interior of the bilhard reads 

j d^g2i^(q3,q2;t')^(q2,qi;t) = j d^g2i^o(q3,q2;t')^o(q2,qi;t) 



+ j d%iro(i?q3,q2;t')^o(i?q2,qi;t) 

± j d%i^o(i?q3,q2;t')^o(q2,qi;t) 

± j d^g2 i^o(q3, q2; t')^o(i?q2, qi; t) . (B.3) 



Using (Ib:2|) and = id turns (iRSj) into 

y d%iro(q3,q2;t')^o(q2,qi;t) + ^ d% iro(q3, i?q2; 0^o(i?q2, qi; t) 

± y d%i^o(i?q3,q2;t')^o(q2,qi;t)± J d% iro(i?q3, i?q2; t')^o(i?q2, qi; 

= J d^g2 i^o(q3, q2; t')^o(q2, qi; t) ± J d^gs i^o(^q3, q2; t')^o(q2, qi; t) . (B.4) 
nuR(n) QuR(n) 

As the intermediate coordinates q2 now are no longer restricted to the interior Q but 
instead can be regarded as running over full space, the convolution property of the free 
propagator can be applied to obtain 

J d^gs i^(q3, q2; t')^(q2, ^ut) = M^is, qi; t + t')± MRqs, qi; t + 1') 

n 

= ir(q3,qi;t + t')- (B.5) 

This shows that the assumption of local flatness of the boundary suffices to give the 
confined propagator the convolution property needed for the convolution formula by 
Weidenmiiller to hold also when using it in combination with the single-particle Weyl 
expansion up to the first boundary correction instead of the exact DOS. 

Appendix C. Calculation via Convolution Formula - confined case 

As in the unconfined case, the many-body DOS can be computed from convolutions of 
single-particle DOS in Weyl expansion. The single-particle Weyl expansion including 
the first boundary correction reads 



i2_ 3 
2 2 



= aE^-^e{E) + bE^-le{E). (C.l) 
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Using (1521) and writing the sum of partitions as ordered tuples yields the sum 



N 



Qconv,i: 



i±i 



N 



1=1 



n 



E n 



Ni,...,Ni=l u}=l 
Y^Nu;=N 



1 \2 



C{E) 



of convolutions C{E) of the form 



C{E) 



dEi . . . dEi 



UJ = 1 



6[e-Y.E. 



Ul=l 



lv=0 ^ 



where 



1 \ o i- 



U]=l 



n 

U]=ly+1 



D _ 3 
2 2 ~ 



with 



/•CO /JX- D \ / ^ D 3\ / ^ 



(C.2) 



(C.3) 



(C.4) 



(C.5) 



The reason for expressing (1C.2P as the sum over ordered tuples is the resulting invariance 
of the sum with respect to relabelling amongst the Ni, . . . ,Ni in the summands. This 
allows to subsume all summands in nL=i Qsp{Eu)/ N^) with certain powers in a and h. 
The number of such equivalent summands is counted by the binomial coefficient in ( 1C.3I) . 
Again the partial sums of energies 



E- J2 



n = / 



(C.6) 



a;=n+l 



are defined. The first ly integrals then are identical to the unconfined case with the 
substitutions I — )■ ly and E ^ ei^, yielding 



X{E) 



r(^ + i) 



r(f)'^ 



dEi^+i...dEi Yl E, 

ti;=iv+l 



i2_ 3 

2 2 



(C.7) 



The /s = I — W remaining integrals are all again of the form (I A. 71) where this time 



D 3 



— |. Also the initial coefficient Ai and exponent si are different, all in all leading 



2 2 



to the following recurrence relation 



_ r(i + r)r(i + .o , _r(^ + i) 

^n+l — T^/fi , ; ^ ' ^1 



" r(2 + r + 

Sn+l = r + Sn + l, 



r(f) 



(C.8) 



A geometrical approach to the mean density of states in many-body quantum systemsAA 



solved by 



Sn+i = Si + n{r + 1) , 

r(l + si) (C.9) 



An+i = AiTil + r) 



r(i + s„+i: 



Recognizing that 



X{E)=Ai,+^E'^s+^e{E), (C.IO) 
and reintroducing the expressions for a and h leads to 



^<-)-Ea)(n^)*-'(iiJ 

tv=0 1^=1 a;=(v+l 



2 2 



(c.ir 

X ys ^^"^^ ^ eiE) 



which, inserted into (1C.2I) . gives exactly the expression (1751) . 

One has to mention that for D = 1 the second term in the single-particle Weyl 
expansion has to be replaced according to 



e{E)^6{goE). (C.12) 



Nevertheless this does not heavily change the obtained result, as in the case ly > the 
integrals involved instead of (]A.7p are of the form 



dx5ix)ic~xr = c^9ic) = c-+-+ Y^^^| + ^^^^ g(c) 



(C.13) 



which leads to a similar recurrence relation. Recognizing that the missing factor of 
r(l + r) in the end is just compensated by the replacement of T{{D — l)/2) in flC.12p 



this shows the validity of all terms with ly > also in the case D = 1. The summand 
corresponding to ly = involves convolutions of Dirac-Delta distributions only, leading 
to the already mentioned replacement rule ( !78l) in the final expression similar to ( 1C.12I) . 

It has also to be mentioned that analogous to the presented calculation it is possible 
to include also higher terms in the single-particle Weyl expansion {D > 2) corresponding 
to corrections originating in the curvature of the boundary. The resulting expressions 
are then involving an additional summation over an index Iq representing the number 
of clusters contributing via curvature corrected propagation. The already time intense 
computation of the coefficients Ci^i^ fl76|) without curvature corrections then becomes 
even worse in the further extended case. Also the property (|83|) had then to be shown 
to link the smooth part of the DOS with the convolution formula of Weidenmiiller. 

As examples (see section [8]) show, the effect of the curvature in singly connected 
billiards generically is not as strong as the effect of the surface correction and might 
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be approximated by simply shifting the energy corresponding to the associated Eq^. 
Therefore and because of the analogous and straightforward but extensive derivation of 
the corresponding DOS, this computation is not explicitly shown here. 
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